#######################
# UNFAIR FIGHTS: POWER ASYMMETRY, NASCENT NUCLEAR CAPABILITY, AND PREVENTIVE CONFLICT 
# SCHUB
# CMPS REPLICATION FILE 
#######################


####################
# CONTENTS OVERVIEW
####################

# (I)  Figure 1 from Paper
# (II) Figures 2-4 from Supporting Information

####################
# (I)  Figure 1 from Paper
####################
library(foreign)
dataraw<-read.dta("ShiftYear.dta")
data<-dataraw[dataraw$forongoing==0,]
dim(data)
head(data)

n<-0.25
cut<-seq(0,1,by=n)
length(cut)
mat<-matrix(NA,ncol=4,nrow=length(cut))
for(i in 1:length(cut)){
#i<-2
mat[i,1]<-cut[i]	
a<-data$mid[data$reldef>cut[i] & data$reldef<=(cut[i]+n)]
mat[i,2]<-mean(na.omit(a))
mat[i,3]<-sd(na.omit(a))
mat[i,4]<-length(na.omit(a))
}
mat

plot((mat[,1]+0.125),mat[,2]*100,xlim=c(0,1),ylim=c(0,19.75),pch=15,xaxt="n",xlab="",yaxt="n",ylab="",cex=1.5)
rug(data$reldef[data$mid==0][seq(1,length(data$reldef[data$mid==0]),by=2)],lwd=0.01,side=1,col="grey60")
rug(data$reldef[data$mid==1],lwd=1.55,side=1,col="black")
abline(v=0,lty=3,lwd=3)
abline(v=0.25,lty=3,lwd=3)
abline(v=0.5,lty=3,lwd=3)
abline(v=.75,lty=3,lwd=3)
abline(v=1,lty=3,lwd=3)
axis(side=1,at=seq(0,1,by=0.25),lab=c("0","25","50","75","100"),mgp=c(0,0.6,0))
mtext("Declining State's Share of Capabilities (%)",outer=T, line=-2.3,side=1,cex=1.2)
axis(side=2,at=seq(0,20,by=4))
mtext("Probability of Conflict (%)",outer=T, line=-1.7,side=2,cex=1.2)
mtext("Descriptive Evidence",outer=T, line=-3.4,side=3,cex=1.3)
mtext("Shifting Power, Power Asymmetry, and Conflict",outer=T, line=-2.0,side=3,cex=1.3)
text(0.36,10.9,"quarter's mean")
segments(.35,(8.99+2.4),.37,(10.39+2.4))
segments(.77,2.45,.77,2,col="black",lwd=1.5)
segments(.77,1.37,.77,0.92,col="grey60",lwd=1.5)
text(0.85,2.26,"= conflict")
text(0.87,1.15,"= no conflict", col="grey40")


####################
# (II)  Figures 2-4 from Supporting Information
####################
library(foreign)
library(polywog)
library(mgcv)

################
###### Figure 2
  
  #### cubic polynomials
dataraw<-read.dta("ShiftYear.dta")  
rjs<-dataraw[dataraw$forongoing==0,]  
dim(rjs)  


rjs$reldef2<-rjs$reldef*rjs$reldef
rjs$reldef3<-rjs$reldef*rjs$reldef*rjs$reldef
def<-glm(mid ~ reldef +reldef2 + reldef3 + cont + pol21 + pol22 + duration3 + pceyrs + pceyrs2 + pceyrs3, data=rjs, family=binomial(link=logit))
summary(def)
test=data.frame(reldef=seq(0,1,by=0.01), reldef2=(seq(0,1,by=0.01)^2), reldef3=(seq(0,1,by=0.01)^3), cont= median(def$model$cont), pol21= 0, pol22= 0, duration3= median(def$model$duration3), pceyrs= median(def$model$pceyrs), pceyrs2= median(def$model$pceyrs2), pceyrs3= median(def$model$pceyrs3))
ft<-predict(def,test,type="response",se.fit=TRUE)
reldef<-seq(0,1,by=0.01)


  
  #### gam
gm<-gam(mid ~ cont + pol21 + pol22 + duration3 + pceyrs + pceyrs2 + pceyrs3 + s(reldef), data=rjs, family=binomial)
summary(gm)
test=data.frame(reldef=seq(0,1,by=0.01), cont= median(gm$model$cont), pol21= 0, pol22= 0, duration3= median(gm$model$duration3), pceyrs= median(gm$model$pceyrs), pceyrs2= median(gm$model$pceyrs2), pceyrs3= median(gm$model$pceyrs3))
fits=predict(gm,newdata=test,type='response',se=T)
predicts=data.frame(test,fits)

  #### PLOT
#pdf("monotonic1.pdf",height=6, width=8)
par(mfrow=c(1,2))
par(oma=c(2.5,3,2,0))
par(mar=c(2,1.5,3,1))
plot(reldef,ft$fit,type="l",lwd=4,ylim=c(0,0.49),xaxt="n",yaxt="n")
points(reldef,(ft$fit-1.96*ft$se.fit),type="l",lwd=3,lty=3,col="grey50")
points(reldef,(ft$fit+1.96*ft$se.fit),type="l",lwd=3,lty=3,col="grey50")
axis(1,at=c(0,0.25,0.5,0.75,1),cex.axis=1.1)
axis(2,at=c(0,0.1,0.2,0.3,0.4,0.5),labels=c(0,10,20,30,40,50),cex.axis=1.1)
mtext("Declining State's Share of Capabilities (%)",outer=T, line=0.6,side=1,cex=1.2)
mtext(cex=1.4,side=2,outer=TRUE,line=1,"Probability of Conflict (%)")
mtext(cex=1.6,side=3,outer=TRUE,line=-0.5,"Assessing Monotonocity (I)")
mtext(cex=1.4,side=3,outer=FALSE,line=0.6,"Polynomials (cubic)")
plot(predicts$reldef,predicts$fit,type="l",ylim=c(0,0.49),lwd=4,yaxt="n",xaxt="n")
points(predicts$reldef,(predicts$fit-1.96*predicts$se.fit),type="l",lwd=3,lty=3, col="grey50")
points(predicts$reldef,(predicts$fit+1.96*predicts$se.fit),type="l",lwd=3,lty=3, col="grey50")
axis(1,at=c(0,0.25,0.5,0.75,1),cex.axis=1.1)
axis(2,at=c(0,0.1,0.2,0.3,0.4,0.5),labels=c(0,10,20,30,40,50),cex.axis=1.1)
mtext(cex=1.4,side=3,outer=FALSE,line=0.6,"Generalized Additive Model")
#dev.off()

################
###### Figure 3
  
  #### binary quarter indicators
# see do file for simulations generating estimates
est<-c(12.1,21.5,25.1,27.7)
lo<-c(5.7,12,14.4,20.1)
hi<-c(21,33,39.3,36.3)
range<-c(.125,.375,.625,.875)

  #### Polywog
pw<-polywog(mid ~ reldef + cont + pol21 + pol22 + duration3 + pceyrs + pceyrs2 + pceyrs3, data=rjs, degree=2, boot=40)
summary(pw)
## borrow the mean estimates for other vars from the gam 
dt=data.frame(reldef=seq(0,1,by=0.01), cont= median(gm$model$cont), pol21=0, pol22=0, duration3= median(gm$model$duration3), pceyrs= median(gm$model$pceyrs), pceyrs2= median(gm$model$pceyrs2), pceyrs3= median(gm$model$pceyrs3))
reldef<-seq(0,1,by=0.01)
fits<-predict(pw,dt,level=0.95,interval=TRUE)




  #### PLOT
#pdf("monotonic2.pdf",height=6, width=8)
par(mfrow=c(1,2))
par(oma=c(2.5,3,2,0))
par(mar=c(2,1.5,3,1))
plot(range,est,ylim=c(0,49),xlim=c(0,1),col="white",yaxt="n",xaxt="n")
segments(range,lo,range,hi,col="grey50",lwd=3)
points(range,est,pch=15,cex=1.5)
axis(1,at=c(0,0.25,0.5,0.75,1),cex.axis=1.1)
axis(2,at=c(0,10,20,30,40,50),labels=c(0,10,20,30,40,50),cex.axis=1.1)
abline(v=0,lty=3,lwd=1)
abline(v=0.25,lty=3,lwd=1)
abline(v=0.5,lty=3,lwd=1)
abline(v=.75,lty=3,lwd=1)
abline(v=1,lty=3,lwd=1)
mtext("Declining State's Share of Capabilities (%)",outer=T, line=0.6,side=1,cex=1.2)
mtext(cex=1.4,side=2,outer=TRUE,line=1,"Probability of Conflict (%)")
mtext(cex=1.6,side=3,outer=TRUE,line=-0.5,"Assessing Monotonocity (II)")
mtext(cex=1.4,side=3,outer=FALSE,line=0.6,"Indicators for Each Quarter")
plot(reldef,fits[,1],ylim=c(0,0.49),type="l",lwd=4,xaxt="n",yaxt="n")
points(reldef,fits[,2],type="l",col="grey50",lty=3,lwd=3)
points(reldef,fits[,3],type="l",col="grey50",lty=3,lwd=3)
mtext(cex=1.4,side=3,outer=FALSE,line=0.6,"Bootstrapped Basis Regression")
axis(1,at=c(0,0.25,0.5,0.75,1),cex.axis=1.1)
axis(2,at=c(0,0.1,0.2,0.3,0.4,0.5),labels=c(0,10,20,30,40,50),cex.axis=1.1)
#dev.off()


################
###### Figure 4
  
  #### cubic polynomials
prd<-read.dta("PRDyear.dta")
nos<-prd[prd$nukeshift==0,]
dim(nos)

reldef<-seq(0,1,by=0.01)
nos$reldef2<-nos$reldef*nos$reldef
nos$reldef3<-nos$reldef*nos$reldef*nos$reldef
def<-glm(mid ~ reldef +reldef2 + reldef3 + cont + pceyrs + pceyrs2 + pceyrs3, data=nos, family=binomial(link=logit))
summary(def)
tes=data.frame(reldef=seq(0,1,by=0.01), reldef2=(seq(0,1,by=0.01)^2), reldef3=(seq(0,1,by=0.01)^3), cont= median(def$model$cont),pceyrs= median(def$model$pceyrs), pceyrs2= median(def$model$pceyrs2), pceyrs3= median(def$model$pceyrs3))
ft<-predict(def,tes,type="response",se.fit=TRUE)

  #### gam [sample from data to get it to run]
snos<-nos[c(seq(1,nrow(nos),by=10)),]
dim(snos)

gm<-gam(force ~  cont + pceyrs + pceyrs2 + pceyrs3 + s(reldef), data=snos, family=binomial)
summary(gm)
test=data.frame(reldef=seq(0,1,by=0.01), cont= 0, pceyrs= 11, pceyrs2= 121, pceyrs3= (121*11))
fits=predict(gm,newdata=test,type='response',se=T)
predicts=data.frame(test,fits)

  #### PLOT
par(mfrow=c(1,2))
par(oma=c(2.5,3,2,0))
par(mar=c(2,1.5,3,1))
plot(reldef,ft$fit,type="l",lwd=4,ylim=c(0,0.005),yaxt="n",xaxt="n")
points(reldef,(ft$fit-1.96*ft$se.fit),type="l",lwd=3,lty=3,col="grey50")
points(reldef,(ft$fit+1.96*ft$se.fit),type="l",lwd=3,lty=3,col="grey50")
axis(1,at=c(0,0.25,0.5,0.75,1),cex.axis=1.1)
axis(2,at=c(0,0.0010,0.002,0.003,0.004,0.005),labels=c(0,0.10,0.20,0.30,0.40,0.50),cex.axis=1.1)
mtext("Declining State's Share of Capabilities (%)",outer=T, line=0.6,side=1,cex=1.2)
mtext(cex=1.4,side=2,outer=TRUE,line=1,"Probability of Conflict (%)")
mtext(cex=1.6,side=3,outer=TRUE,line=-0.5,"Relationship without Shifting Power")
mtext(cex=1.4,side=3,outer=FALSE,line=0.6,"Polynomials (cubic)")
plot(predicts$reldef,predicts$fit,type="l",ylim=c(0,0.005),lwd=4, yaxt="n",xaxt="n")
points(predicts$reldef,(predicts$fit-1.96*predicts$se.fit),type="l",lwd=2,lty=3)
points(predicts$reldef,(predicts$fit+1.96*predicts$se.fit),type="l",lwd=2,lty=3)
axis(1,at=c(0,0.25,0.5,0.75,1),cex.axis=1.1)
axis(2,at=c(0,0.0010,0.002,0.003,0.004,0.005),labels=c(0,0.10,0.20,0.30,0.40,0.50),cex.axis=1.1)
#rug(rjs$reldef,lwd=0.01)
mtext(cex=1.4,side=3,outer=FALSE,line=0.6,"Generalized Additive Model")  


